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Abstract 

An  approach  to  modeling  the  impact  of  disturbances  in  an  agricultural  production 
network  is  presented.  A  stochastic  model  and  its  approximate  deterministic  model  for 
averages  over  sample  paths  of  the  stochastic  system  are  developed.  Simulations,  sensi¬ 
tivity  and  generalized  sensitivity  analyses  are  given.  Finally,  it  is  shown  how  diseases 
may  be  introduced  into  the  network  and  corresponding  simulations  are  discussed. 
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1  Introduction 

The  current  production  methods  for  livestock  follow  the  just-in-time  philosophy  of  man¬ 
ufacturing  industries.  Feedstock  and  animals  are  grown  in  different  areas.  Animals  are 
moved  from  one  farm  to  another,  depending  on  their  age.  Unfortunately,  shocks  propagate 
rapidly  through  such  systems;  an  interruption  to  the  feed  supply  has  a  much  larger  impact 
when  farms  have  minimal  surplus  supplies  in-stock  than  when  they  have  large  reserves.  The 
just-in-time  movement  of  animals  between  farms  serves  as  another  vulnerability.  Stopping 
movement  of  animals  to  and  from  a  farm  with  animals  infected  by  a  disease  will  have  effects 
that  quickly  spread  through  the  system.  Nurseries  supplying  the  farm  will  have  nowhere  to 
send  their  animals  as  they  grow  up.  Finishers  and  slaughterhouses  supplied  by  the  farm  will 
have  their  supply  interrupted. 

The  devastating  Foot  and  Mouth  Disease  (FMD)  that  hit  the  United  Kingdom  (UK)  in 
2001  lead  to  the  slaughter  of  millions  of  animals.  The  outbreak  shook  many  western  nations 
as  they  watched  a  nation  with  an  advanced  animal  health  surveillance  and  response  system 
fail  to  get  FMD  under  control,  in  part  because  they  were  unable  to  mount  a  rapid  response  in 
the  face  of  modern  agricultural  marketing  systems  [15].  In  an  effort  to  eradicate  the  disease, 
the  marketing  channels  were  stopped,  putting  uninfected  producers  in  the  situation  of  having 
no  income  to  maintain  their  livestock  and  no  means  to  move  them  to  locations  where  feed, 
shelter,  and  other  support  were  available.  As  a  result,  between  six  million  and  ten  million 
animals  were  destroyed  in  the  UK  over  seven  months,  with  over  one-third  of  those  animals 
being  destroyed  for  welfare  reasons  [41],  Two  years  after  the  outbreak,  animal  agriculture 
in  the  UK  was  still  declining,  a  chilling  postscript  to  the  widespread  infrastructure  damage 
FMD  had  wrought  on  the  nation  [37]. 

More  recently,  the  world  has  witnessed  the  apparent  failure  of  widespread  national  and 
international  plans  for  using  animal  destruction  to  stem  the  spread  of  the  highly  pathogenic 
H5N1  strain  of  avian  influenza.  In  a  process  frighteningly  reminiscent  of  the  UK  FMD 
experience,  the  programs  have  also  allowed  domestic  markets  within  and  beyond  affected 
countries  to  suffer.  Global  consumption  of  poultry  has  dropped  sufficiently  to  cause  US 
domestic  producers  (e.g.,  Tyson,  Pilgrims  Pride,  et  ah)  to  absorb  decreased  demand  and 
decreased  prices.  This  drop  has  translated  to  decreases,  as  well,  in  non-poultry  markets, 
exacerbating  the  market  effects  of  a  disease  not  even  present  in  the  western  hemisphere 
[19].  It  has  become  painfully  apparent  that  in  the  large-scale,  interdependent,  and  highly 
mobile  animal  agriculture  industry  of  the  USA,  the  unintended  consequences  and  market 
ripple  effects  of  a  disease  incursion  into  our  system  could  be  even  more  severe  than  what  was 
witnessed  in  the  UK  in  2001  and  across  Europe  in  2005-6,  and  could  induce  decision  makers 
to  call  for  even  more  draconian  measures  than  previously  seen.  What  is  needed  is  a  new  view 
of  how  our  emergency  response  programs  might  affect  modern  animal  agriculture,  a  view 
that  allows  workers  to  assess  the  potential  for  other  prevention  strategies  and  responses.  The 
view  should  also  allow  analysts  to  identify  bottlenecks  in  the  food  and  feed  supply  chain, 
and  to  test  potential  mitigation  tools,  procedures,  and/or  practices  to  increase  the  resilience 
of  animal  agriculture  to  catastrophic  animal  diseases. 
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This  paper  presents  initial  statistical  and  mathematical  modeling  ideas  to  address  the 
above  issues,  using  the  North  Carolina  swine  industry’s  potential  response  to  FMD  as  an 
example.  We  focused  our  attention  on  the  North  Carolina  swine  industry  because  it  is  both 
the  second  largest  swine  industry  in  the  United  States,  and  is  local  to  us.  Our  goal  was  to 
develop  a  model  that  could  be  used  to  investigate  how  small  perturbations  to  the  agricultural 
supply  system  would  affect  its  overall  performance.  A  hurricane  that  throttles  inter-farm 
transportation  for  a  short  duration,  or  a  disease  outbreak  that  spreads  through  distribution 
channels  are  example  causes  of  the  perturbations  of  interest.  In  the  former  case,  the  just-in- 
time  delivery  systems  may  not  provide  enough  slack  to  absorb  the  shock.  In  the  latter  case, 
strategies  that  involve  destruction  of  all  livestock  in  an  infected  branch  of  the  system  may 
be  overly  harsh;  a  more  moderate  response  may  be  as  effective  without  the  high  toll  on  the 
infrastructure. 

We  model  a  simplified  swine  production  network  in  North  Carolina  containing  four  levels 
of  production  nodes:  growers/sows  ( Node  1),  nurseries  ( Node  2),  finishers  ( Node  3),  and 
processing  plants/slaughter  houses  ( Node  4).  At  grower  or  sow  farms  ( Node  1),  the  new 
piglets  are  born  and  typically  weaned  three  weeks  after  birth.  The  three-week  old  piglets  are 
then  moved  to  the  nursery  farms  ( Node  2)  to  mature  for  another  seven  weeks.  They  are  then 
transferred  to  the  finisher  farms  ( Node  3),  where  they  grow  to  full  market  size,  which  takes 
about  twenty  weeks.  Once  they  reach  market  weight,  the  matured  pigs  are  moved  to  the 
processors  (slaughter  houses)  ( Node  4).  Pork  products  then  continue  through  wholesalers  to 
consumers.  There  are  also  several  inputs  to  the  system  which  we  will  not  consider,  such  as 
food,  typically  corn  grown  in  the  midwest.  There  are  several  types  of  breeding  farms  where 
pure-bred  stock  are  raised;  these  are  typically  crossed  to  produce  hybrid  strains  for  pork 
production. 

Our  paper  is  organized  as  follows.  In  Section  2,  we  formulate  a  nonlinear  stochastic 
model  for  our  agricultural  network  and  show  how  it  can  be  converted  to  an  equivalent  (in 
a  sense  made  precise  below)  deterministic  differential  equation  model.  This  deterministic 
model  readily  lends  itself  to  simulations  and  sensitivity  analysis  techniques.  In  Section  3 
we  present  numerical  simulations  of  the  production  model  (without  perturbations  such  as 
infectious  disease),  and  carry  out  a  sensitivity  analysis  of  the  model.  Simulations  of  the  model 
in  the  presence  of  an  infectious  disease  are  presented  in  Section  4.  Finally,  in  Section  5  we 
give  our  conclusions  and  remarks  for  future  work. 

In  addition  to  the  development  of  models  for  a  typical  production  network  permitting 
perturbations,  a  significant  contribution  in  this  paper  is  the  demonstration  of  stochastic, 
mathematical  and  computational  methodology  that  is  available  to  domain  scientists,  statisti¬ 
cians  and  applied  mathematicians  working  in  a  concerted  team  effort  on  complex  problems 
of  the  type  exemplified  here.  The  co-authors  of  this  paper  constituted  such  a  team  orga¬ 
nized  under  the  auspices  and  with  the  support  of  the  Statistical  and  Applied  Mathematical 
Sciences  Institute  (SAMSI)  as  a  year  long  working  group  in  its  recent  research  program  on 
National  Defense  and  Homeland  Security. 
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2  Modeling 

We  consider  stochastic  models  to  track  an  agricultural  network.  We  are  interested  in  how  the 
parameters  used  in  the  model  affect  the  overall  capacity  of  a  network,  and  how  one  discerns 
the  existence  and  location  of  any  bottlenecks.  With  deterministic  models,  one  can  answer 
the  first  question  with  a  sensitivity  analysis.  Thus,  after  developing  a  typical  stochastic 
production  model,  we  also  show  how  to  obtain  its  deterministic  approximation.  We  then 
demonstrate  how  to  superimpose  a  simple  contagious  disease  model  on  the  production  model 
that  allows  simulation  of  dynamics  and  spread  of  FMD  through  a  production  chain. 

2.1  Basic  Model 

We  consider  a  simplified  swine  production  network  with  four  aggregated  nodes:  sows  ( Node  1), 
nurseries  ( Node  2),  finishers  ( Node  3),  and  slaughter  houses  ( Node  4).  Our  goal  is  to  study 
the  effects  of  perturbations  within  the  network.  This  can  be  done  either  by  affecting  the 
nodes  or  the  transitions  between  nodes  directly  or  indirectly.  For  instance,  a  problem  with 
the  breeding  farms  would  result  in  a  reduction  of  sows  available  for  producing  new  piglets. 
This  would  result  in  a  reduced  rate  of  transition  from  Node  1  to  Node  2,  since  we  could  not 
grow  as  many  piglets.  We  could  then  track  the  effect  of  this  through  our  network. 

Although  unavoidable  in  actual  production  processes,  we  assume  in  our  example  that 
there  are  no  net  losses  in  the  network  (i.e.,  the  total  number  of  pigs  in  the  network  remains 
constant)  and  that  the  only  deaths  occur  at  the  slaughter  houses.  Thus  we  assume  that  the 
number  of  processed  pigs  per  day  at  the  slaughter  houses  is  equal  to  the  number  of  newborn 
piglets  per  day  at  the  growers.  We  can  model  reduced  birth-rates  by  reducing  the  rate  at 
which  piglets  move  to  the  nurseries.  This  leads  us  to  deal  with  a  closed  network.  We  note 
that  this  approximation  is  realistic  when  the  network  is  efficient  and  operates  at  or  near  full 
capacity  (i.e,  when  the  number  of  animals  removed  from  the  chain  are  immediately  replaced 
by  new  production/growth,  avoiding  significant  idle  times).  Our  closed  network  model  for 
the  swine  production  is  summarized  schematically  in  Figure  1. 

Sows  Nursery  Finisher  Slaughter 


L-*-  N 1  )  - ►  y  N2  )  - ►  N3  ; - ►  \  N4  ) 

Figure  1:  Aggregated  agricultural  network  model. 

Each  node  with  corresponding  population  number  Nt,  i  —  1, . . . ,  4,  in  Figure  1  represents 
an  aggregation  of  all  the  production  units  corresponding  to  that  level  in  the  production 
network.  Given  a  specific  production  network,  any  of  the  four  levels  of  the  chain  may  be 
broken  into  its  constituent  units  (e.g.,  farms),  and  analyzed  in  detail  as  a  separate  sub¬ 
network.  The  directed  edges  between  the  nodes  represent  the  movement  of  the  pigs  through 
the  network.  The  rate  is  determined  by  the  pigs’  residence  time,  the  number  of  pigs  at  each 
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node,  and  the  capacity  constraints  at  the  corresponding  nodes.  Let  L*  denote  the  capacity 
constraint  at  node  i,  i  =  1, ...  ,4.  Since  we  have  a  closed  network,  it  is  assumed  that  there 
is  no  capacity  constraint  at  Node  1  and  therefore  we  take  L4  =  oo.  We  also  define  Sm  to 
be  the  maximum  exit  rate  at  Node  4,  i.e.,  the  maximum  killing  capacity  at  the  slaughter 
house.  The  residence  times  at  each  node,  together  with  the  capacity  constraints  and  the 
slaughterhouse  killing  capacity,  based  on  very  rough  estimates  of  swine  production  in  North 
Carolina  [1],  are  given  in  Table  1. 


Name 

Sows 

Nursery 

Finisher 

Slaughter 

Node 

N1 

N2 

N3 

N4 

Piglet  residence  time  (days) 

21(jVl-UV2) 

4:9(W2— >iV3) 

140(tv3^7V4) 

l(jV4^AU) 

Assumed  capacity  (in  thousands) 

00 

825 

2,300 

20 

Table  1:  Network  parameters  based  on  swine  production  in  NC. 


2.2  Stochastic  and  Deterministic  Models 

We  model  the  evolution  of  the  food  production  network  shown  in  Figure  1  as  a  continuous 
time  discrete  state  density  dependent  jump  Markov  Chain  (MC)  [3,  21]  with  a  discrete  state 
space  embedded  in  an  M4  non- negative  integer  lattice  £.  The  state  of  this  MC  at  time  t  is 
denoted  by  X(f)  =  (Xi(f), . . . ,X4 (t) ) ,  where  Xi(t)  is  the  number  of  pigs  at  node  i  at  time 
t,i  =  1, ...  ,4. 

The  rates  of  transition  of  X(t)  are  nonlinear  functions  A*  :  C  — >  [0,  00)  for  i  —  1, . . . ,  4, 
and  for  xgf  are  given  by: 


Ai(x) 

:=  gi(o:i  -  l,x2  +  l,x3,  o;4) 

=  k1x1(L2-  x2)  + 

A2(x) 

:=  q2(xi,x2  ~  1,0:3  +  1,0:4) 

=  k2x2(L3-x3)  + 

*3(x) 

:=  qz{x  1,0:2, 0:3  -  1,0:4  +  1) 

=  k3x3(L4-  x4)+ 

A4(x) 

:=  g4(o:i  +  1,  x2,  x3,  x4  -  1) 

=  k4  min(o:4,  Sm) 

(1) 

where  kt,i  =  1,...,4,  is  proportional  to  the  service  rate  at  node  i ;  =  2,3,4,  is  the 

buffer  size  (capacity  constraint)  at  node  i  and  Sm  is  the  slaughter  capacity  at  node  4  as 
discussed  above.  For  any  real  z,  the  symbol  (z)+  is  defined  as  the  non-negative  part  of  z, 
i.e.,  (z)+  =  max(z,  0).  Then  qi(x\  —  l,x2  +  1,0:3,  £4)  is  given  by 

Pr[X(f  +  h)  =  (xi  —  1,0:2  +  l,:r3,  z4)|X(f)  =  (x1,x2,  x3,x4)} 

qi{x1-l,x2  +  l,x3,x4)  =  Inn  - - - 

h — >0+  h 


The  other  g*  are  given  similarly. 
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Let  Ri(t)  i  —  1, . . . ,  4,  denote  the  number  of  times  that  the  i ’th  transition  occurs  by  time 
t.  Then  Ri  is  a  counting  process  with  intensity  A *(X(t)),  and  the  corresponding  stochastic 
process  can  be  dehned  by 

Ri(t)  =  Yi([t  A,(X(S))cfe),  i  =  (2) 

Jo 

where  the  Yi  are  independent  unit  Poisson  processes.  That  is,  sample  paths  ?y(f)  of  Ri(t) 
are  given  in  terms  of  sample  paths  x(f)  of  X(t)  by 

n(t)  =  Yi((  A  j(x(s))ds),  i  =  1 . . . ,  4.  (3) 

Jo 

We  write  Ri  in  this  form  to  illustrate  that  \  is  a  rate  of  the  corresponding  counting  process. 

Let  ej,  i  —  1, . . . ,  4,  be  standard  basis  vectors  of  M4  and  define,  for  i  incremented  by  one 
modulo  4,  the  vectors 

Vi  ® (moci4)  *  1)2,...,  4, 

which  represent  the  vector  of  changes  in  system  counts  at  ?’th  transition.  We  write  the  state 
of  the  system  at  time  t  as 

X(t)  =  X(0)  +  ^  Ri(t)ui  =  X(0)  +  uR(t),  (4) 

i 

where  v  is  the  matrix  with  rows  given  by  the  vL1  and  R (f)  is  the  (column)  vector  with 
components  Ri(t).  In  the  chemical  literature,  the  matrix  v7  is  often  referred  to  as  the 
stoichiometric  matrix  [29].  More  specifically,  we  have 

X1(t)  =  X^o)  -  R^t)  +  R4(t) 
x2(t)  =  X2(0)  +  Ri(t)  -  Rv(t) 

X3(t)  =  X3(0 )  +  R2(t)  -  R3(t) 

X4(t)  =  X4(0)+R3(t)-R4(t).  (5) 

The  above  system  typically  cannot  be  solved  for  a  stationary  distribution  and  an  empirical 
approach  based  on  the  so-called  Gillespie  algorithm  [29]  can  be  used  to  investigate  the  long 
term  behavior  of  the  system  (see  Section  3.2).  The  approximate  large  population  behavior 
of  an  appropriately  scaled  system  may  be  also  analyzed  via  macroscopic  deterministic  rate 
equations  as  we  shall  explain  next  (the  original  theory  is  due  to  Kurtz  and  is  discussed  in 
[21]  and  the  references  therein). 

Let  N  be  the  total  network  or  population  size.  If  N  is  known  we  may  consider  the 

animal  units  per  system  size  or  the  units  concentration  in  the  stochastic  process  C N  (t)  = 

X(t)/N  with  sample  paths  cN(t).  For  large  systems  this  approach  leads  to  a  deterministic 
approximation  (obtained  as  solutions  to  the  system  rate  equation  dehned  below)  to  the 
stochastic  equation  (4),  in  terms  of  c(f),  the  large  sample  size  average  over  sample  paths  or 
trajectories  c N(t)  of  C N(t). 
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We  rescale  the  rate  constants  kt,  Lt  and  Sm  as  follows: 

^4  =  ki,  Ki  =  Nki,  i  =  1, 2  or  3, 

sm  =  Sm/N,  u  =  U/N.  (6) 

According  to  Equation  (1),  this  rescaling  implies  that 

Aj(x)  =  m  Xi(Li+1  -  Xi)+/N  =  N Ki  cf  (li+1  -  cf )+  i  =  1,  2, 3, 

and 

A4 (x)  =  k4  min(x4,  Sm)  =  Nk^  min(cf ,  sm). 

Recall  that  for  large  N  the  Strong  Law  of  Large  Numbers  (SLLN)  for  the  Poisson  Process 
Y  implies  Y(Nu)/N  &  u  [30].  One  can  use  this  fact,  along  with  the  rescaling  of  the  constants 
as  given  above,  to  argue  that  sample  paths  rft)  for  the  counting  process  (2)  defined  in  terms 
of  the  sample  paths  x(f)  or  c N(t)  =  x(t)/AT  may  be  approximated  for  large  N  in  terms  of  the 
deterministic  variables  c(f),  the  averages  over  sample  paths  or  trajectories  cN(t )  of  C N(t), 
by 


„(A0 


(<)  =  fri (<)  =  fY'(J  Ai(x(s))*) 

=  fY<(N  /AiCf(S)(/j+I  -  c'X1(S))+  ds) 

~  /  KiCi(s)(li+1  —  Ci+i(s))+ds  fori  =  1,2, 3, 


(7) 


and  similarly 


r[N\t)  =  jj ri(t)  ~  min(c4(s),  sm)  ds. 


For  a  full  and  rigorous  discussion  of  this  “approximation  in  mean”,  see  Chapters  6.4  and  11 
of  [21]  and  Chapter  5  of  [3].  The  averages  c (t)  satisfy  a  system  of  deterministic  ordinary 
differential  equations  which  can  be  heuristically  derived  by  beginning  with  Equation  (5). 
Upon  dividing  both  sides  of  each  equation  by  N  and  applying  the  above,  we  obtain  the  rate 
equations,  i.e. ,  the  system  of  integral  equations  approximating  via  the  SLLN  the  original 


stochastic  system,  as  follows: 


c?(t) 

c?(t) 

c?(t) 


ci(0)  -  r[N\t)  +  r(4N\t ) 

ci(0)  -  /  KiCi(s)(Z2  -  c2(s))+ds  +  K4min(c4(t),sm) 

Jo 

c2  (o)  +  r[N)  (t)  -  r^JV)  (£) 

c2(0)  -  /  K2C2(s)(/3  -  C3(s))+ Vs  +  /  KiCi(s)(Z2  -  c2(s))+Vs 
Vo  Vo 

c3(0)  +  rfV)(£)  -rf°(f) 

c3(0)  -  /  k3c3(s)(/4  -  c4(s))+  Vs  +  /  k2c2(s)(/3  -  c3(s))+ Vs 
Vo  Vo 

c4(0)  +  r^f)  - 

c4(0)  +  /  re3c3(s)(Z4  -  c4(s))+Vs- re4min(c4(t),sm). 


(8) 


./o 

Upon  approximating  the  cfr(t)  on  the  left  above  by  the  Ci(t )  and  differentiating  the 
resulting  equations,  we  find  that  the  integral  equation  system  is  equivalent  to  a  system  of 
ordinary  differential  equations  for  c(t)  G  M4  given  by 


KiCi(t)(l2  -  c2(f))+  +  rc4min(c4(f),  sm) 

«2C2(t)(73  -  c3(f))+  +  «iCi(t)(Z2  -  c2(t))+ 

«3C3(t)(/4  -  C4(t))  +  +  K2C2(t)(l3  -  C3(t))  + 

Ac4min(c4(f), sm)  +  K3c3(t)(Z4  -  c4(t))+  (9) 

with  the  initial  conditions  c(0)  =  Co-  As  we  shall  see  in  the  next  section,  solutions  of  these 
equations  yield  quite  good  approximations  to  the  sample  paths  of  the  stochastic  system. 

We  remark  that  the  product  nonlinearities  Xi  (Uj+i  —  £i+i)+  of  (1)  where  transportation 
occurs  more  rapidly  the  further  the  node  level  is  from  capacity  (i.e. ,  the  system  reacts  more 
rapidly  to  larger  perturbations  from  capacity)  are  only  one  possible  form  for  these  terms. 
One  could  also  reasonably  argue  for  alternative  terms  of  the  form  Xi  Xi+i  where  Xi+ 1  is  the 
characteristic  function  for  the  set  {(Lj+ 1  —  Xi+\)  >  0}  so  that  the  transportation  rate  from 
a  node  depends  only  on  the  number  available  at  that  node  so  long  as  capacity  at  the  next 
node  has  not  been  reached.  We  remark  that  in  this  case  the  sensitivity  analyses  below  are 
more  difficult  due  to  a  lack  of  continuity  of  the  dynamics  in  the  system  equations. 


Vci(f) 

dt 

dc2(t ) 
dt 

dc3(t ) 
dt 

Vc4(t) 

dt 
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3  Computations  and  Model  Comparison 

3.1  Model  Parameter  Values 


In  order  to  carry  out  numerical  simulations  and  to  compare  the  results  of  the  stochastic  and 
deterministic  models  (equations  (5)  and  (9),  respectively),  we  must  choose  reasonable  values 
for  all  model  parameters.  We  note  that  our  paper  focuses  on  methodological  issues  and,  for 
confidentiality  and  proprietary  reasons,  only  limited  information  on  the  swine  production 
network  was  available  to  us.  Thus  some  of  these  parameter  values  may  be  only  rough 
approximations  to  those  that  might  be  obtained  using  inverse  problem  techniques  with  data 
from  actual  production  networks  [1].  Consequently,  the  subsequent  discussions  in  this  paper 
are  in  no  way  an  attempt  to  validate  the  above  models.  Nonetheless,  we  believe  that  the 
order-of-magnitude  approximate  parameter  values  we  are  using  here  are  sufficient  to  allow 
us  to  develop  and  demonstrate  effective  use  of  methods  and  techniques  which  could  be  used 
with  actual  production  network  based  parameters. 

The  parameters  of  the  stochastic  model,  with  the  exception  of  the  transition  rate  con¬ 
stants  ki ,  are  given  in  Table  1.  From  the  expressions  for  the  transition  rates  (1),  we  see  that 
the  residence  times,  t;,  that  pigs  spend  at  node  i  are  given  by 


U 


1  1 

ki  A.j_|_i(f)) 


for  i 


1,  2  or  3, 


(10) 


and  fq  =  l/k±  =  1.  As  discussed  above,  the  nonlinear  form  of  the  transition  rates  (1)  means 
that  the  residence  time  at  a  given  node  depends  on  how  far  the  following  node  is  below  its 
capacity.  Consequently,  we  determine  the  kt  by  assuming  that  the  given  residence  times 
pertain  to  the  network  in  its  equilibrium  state. 

Considering  the  deterministic  model  equations  (9),  we  see  that,  if  there  is  to  be  a  flow 
through  the  system,  the  equilibrium  population  sizes  N*  =  Nc*  at  each  node  must  be  less 
than  the  capacities  of  the  nodes.  It  is  then  straightforward  to  see  that  the  equilibrium 
numbers  of  individuals  at  each  of  the  first  three  nodes  are  proportional  to  the  This 
makes  intuitive  sense  since  there  is  no  loss  as  individuals  move  between  nodes,  and  so,  at 
equilibrium,  the  relative  residence  times  must  equal  the  relative  numbers  of  individuals  at 
the  nodes.  This  argument  need  not  apply  to  the  slaughter  node,  however,  since  individuals 
will  spend  longer  there  than  the  specified  one  day  residence  time  if  the  equilibrium  value  N4 
is  greater  than  Sm.  The  flow  rate  from  node  four  back  to  node  one  is  the  smaller  of  7V|  and 
Sm  and  so  we  have  that 


(n*,n;,n*,n*4 


{tiNl  t2Nl  t3N*4,N*4)  if  N4  <  Sm 
(tiSm,t2Sm,t3Sm,N4)  otherwise. 


(11) 


Notice  that  solving  for  the  equilibrium  of  the  deterministic  model  does  not  give  us  the  value 
of  N4:  since  the  network  is  closed,  the  total  size  of  the  population  is  equal  to  its  value  at 
the  initial  time.  The  values  of  the  ki,  for  %  —  1,  2  and  3,  are  then  given  by 


U  (Li+i  -  JV*+1) ' 


(12) 
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The  parameter  values  and  the  initial  states  for  the  system  (5)  are  tabulated  in  Table  2. 


Parameters 

Definition 

Values 

Units 

k\ 

transition  rate  at  node  1 

1/(90-21) 

1/ (pigs-days) 

k2 

transition  rate  at  node  2 

1/(200-49) 

1/ (pigs-days) 

h 

transition  rate  at  node  3 

1/(5  • 140) 

1/ (pigs-days) 

k4 

transition  rate  at  node  4 

1 

1 / days 

u 

capacity  node  1 

00 

pigs 

L2 

capacity  node  2 

825 

pigs 

L3 

capacity  node  3 

2300 

pigs 

u 

capacity  node  4 

20 

pigs 

Sm 

slaughter  capacity 

440 

pigs 

XM 

initial  condition 

800 

pigs 

X2(t0) 

initial  condition 

700 

pigs 

X3(t0) 

initial  condition 

1500 

pigs 

x4(t0) 

initial  condition 

165 

pigs 

Table  2:  Aggregated  agricultural  network  model:  Parameters  for  stochastic  simulations, 
together  with  our  chosen  initial  conditions.  All  numbers  of  pigs  are  given  in  thousands  here. 

To  obtain  the  parameters  we  use  in  our  deterministic  simulations  we  simply  rescale  the 
parameters  in  Table  2  by  the  total  network  size  N  =  Xi(t0)  +  X2(to)  +  A i3(t0)  +  A"4(t0)  = 
3, 165,  000,  using  equation  (6)  and  Cj(f0)  =  Xi(t0)/N  for  i  —  1, . . . ,  4.  The  results  are  given 
below  in  Table  3. 
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Parameters 

Definition 

Values 

Units 

scaled  rate  at  node  1 

1.674 

1 / days 

^2 

scaled  rate  at  node  2 

0.323 

1 / days 

K  3 

scaled  rate  at  node  3 

4.521 

1 / days 

K4 

scaled  rate  at  node  4 

1 

1 / days 

k 

scaled  capacity  at  node  2 

2.607-  lO”1 

dimensionless 

h 

scaled  capacity  at  node  3 

7.267-  10"1 

dimensionless 

u 

scaled  capacity  at  node  4 

6.3  •  10~3 

dimensionless 

$171 

scaled  slaughter  capacity 

1.390  ■  10"1 

dimensionless 

ci(0) 

scaled  initial  condition 

2.528  ■  10"1 

dimensionless 

c2(0) 

scaled  initial  condition 

2.212  •  10”1 

dimensionless 

c3(0) 

scaled  initial  condition 

4.739  •  lO”1 

dimensionless 

c4(0) 

scaled  initial  condition 

5.21  •  10"2 

dimensionless 

Table  3:  Aggregated  agricultural  network  model:  Rescaled  parameter  values  and  initial 
conditions  for  the  deterministic  model. 
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3.2  Stochastic  Simulations 

The  standard  method  for  the  stochastic  simulation  of  the  discrete  state  continuous  time 
Markov  Chain  of  the  type  considered  here  is  based  on  a  standard  Monte  Carlo  algorithm, 
also  known  as  the  Gillespie  algorithm  [29].  This  algorithm  is  described  below: 

1.  For  a  given  state  of  the  system  x,  compute  Aj(x)  for  i  —  1 . . . ,  M  (in  our  case  M  =  4). 

2.  Calculate  the  summation  of  the  rates  A  =  YliLi  Aj(x)  and  simulate  the  time  until  the 
next  transition  by  drawing  from  an  exponential  distribution  with  mean  1/A. 

3.  Simulate  the  transition  type  IZx  €  {1, ...  ,4}  by  drawing  from  the  discrete  distribution 
with  P(JZx  =  i)  =  Aj(x)/A. 

4.  Update  the  system  state  x  and  repeat. 

Using  the  above  algorithm  implemented  in  the  statistical  software  R  [40] ,  we  carried  out 
numerous  simulations  for  the  model  (5)  with  the  initial  conditions  and  values  for  parameters 
q*  —  (&!■■■  j  &4,  Sm,  L2,  L3,  L4)  given  in  Table  2. 
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Figure  2:  Simulation  of  the  stochastic  model  (5)  of  the  food  production  chain  via  the  Gillespie 
algorithm.  Parameter  values  are  as  given  in  Table  2;  N=3, 165,000. 

A  sample  of  five  realizations  is  plotted  in  Figure  2.  Note  that  the  realizations  exhibit 
very  little  visible  differences.  However,  when  one  carries  out  the  simulations  for  a  smaller 
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system  ( N  =  3, 165  pigs  with  the  parameters  in  Table  2  scaled  accordingly),  the  variations  are 
readily  visible  as  can  be  seen  in  Figure  3.  We  also  remark  that  these  two  figures  offer  graphic 
depictions  of  the  approximation  theory  discussed  in  Section  2.2  where  in  the  case  of  very 
large  N  one  can  cannot  distinguish  between  the  stochastic  simulations  and  the  corresponding 
deterministic  simulations  for  the  sample  path  averages. 
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Figure  3:  Simulation  of  the  stochastic  model  (5)  with  N=3,165  and  parameter  values  from 
Table  2  scaled  accordingly. 


Stochastic  Simulations 


3.3  Deterministic  Simulation 

We  numerically  solved  the  ODE  system  (9)  using  the  odel5s  solver  in  Matlab.  We  fixed 
the  parameters  q  =  q*  to  be  the  same  as  in  the  stochastic  simulation  above  (but  scaled 
as  in  Table  3)  and  graph  the  solution  of  the  rate  equations  (9)  for  t  E  [0, 100]  in  Figure  4. 
In  Figure  4(left),  we  plot  the  numerical  solution  of  the  concentrations  in  system  (9).  To 
facilitate  comparison  with  the  MC  realizations  plotted  in  Figure  2,  we  also  depict  the  rescaled 
quantities  Nt(t)  =  Nc^t),  which  provide  approximations  to  averages  over  sample  paths 
of  Xi(t)  in  Figure  4(right).  As  expected,  we  find  that  the  stochastic  and  deterministic 
computations  provide  similar  numerical  results  with  the  realizations  fluctuating  about  the 
solution  of  the  deterministic  system  for  the  averages  as  predicted  by  the  theory. 

In  order  to  observe  the  finer  dynamics  in  the  network,  we  plot  the  solution  of  the  network 
on  a  smaller  time  scale  in  Figure  5 (left) .  We  find  that  our  model  solutions  approach  the 
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steady  states  rather  quickly.  All  components  remain  stable  thereafter  suggesting  that  the 
steady  states  are  (at  least  locally)  asymptotically  stable  (this  can  be  verified  with  analytical 
arguments).  We  believe  that  this  behavior  of  the  model  describes  the  food  production 
network  realistically  when  it  is  uninterrupted  by  external  events.  In  Figures  5(right),  6(left) 
and  6(right)  we  observe  similar  behavior  in  the  food  production  chain  for  different  values 
of  the  parameters  L4  and  Sm.  As  one  can  observe  in  these  figures,  when  the  value  of  Sm 
is  sufficiently  large,  the  state  N4l  which  is  related  to  the  replenishment  of  the  network, 
will  never  reach  Sm,  making  the  slaughter  capacity  constraint  inactive  in  the  production 
system.  Only  when  Sm  is  smaller  than  a  certain  critical  value  will  it  be  active  and  in  this 
case  we  observe  accumulation  of  animals  in  the  slaughter  house  (e.g.,  see  Figure  6(right)). 
These  calculations  along  with  numerous  others  we  carried  out  suggest  reasonable  stability 
properties  of  the  production  chain  in  the  absence  of  any  interventions  such  as  FMD  (we  will 
investigate  such  disturbances  below). 


Ag  Concentration  System  Ag  ODE  Solution 


Figure  4:  Numerical  solution  of  the  deterministic  (rate  equation)  system  (9)  for  t  G  [0, 100]. 
Parameter  values  are  as  given  in  Table  3. 
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Ag  ODE  Solution  Ag  ODE  Solution  (Sm=440,L4=1 000) 


Figure  5:  (left)  Numerical  solution  of  the  corresponding  deterministic  system  for  t  E  [0,50]. 
(right)  Numerical  solution  of  the  deterministic  system  (9)  for  t  E  [0,  210]  with  L4  =  1000 
and  Sm  =  440.  All  other  parameter  values  are  as  given  in  Table  3. 


Ag  ODE  Solution  (Sm=1 00,L4=1 000) 


time  (days) 


Figure  6:  (left)  Numerical  solution  of  the  deterministic  system  (9)  for  t  E  [0, 210]  with 
L4  =  1000  and  Sm  =  100.  (right)  Numerical  solution  of  the  deterministic  system  (9)  for 
t  E  [0,  210]  with  L4  =  1000  and  Sm  =  5.  All  other  parameter  values  are  as  given  in  Table  3. 


3.4  Sensitivity  Analysis 

In  this  section,  we  perform  a  sensitivity  analysis  of  the  deterministic  model  (9),  investigating 
how  much  the  solution  of  the  system  changes  when  the  rates  «*,  the  capacities  /*,  or  the  initial 
conditions  Coj,  %  —  1, . . . ,  4  change.  This  analysis  will  be  used  to  identify  the  parameters  and 
the  initial  conditions  to  which  the  system  is  the  most  and  least  sensitive. 

A  second  issue  we  address  here,  which  is  of  great  interest  for  inverse  or  parameter  esti- 
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mation  problems  in  a  typical  nonlinear  regression  model,  is  the  sensitivity  of  the  parameter 
estimates  with  respect  to  the  data  measurements.  We  carry  out  this  analysis  by  means  of 
the  generalized  sensitivity  functions  (GSF)  recently  introduced  by  Thomaseth  and  Cobelli 
[44];  these  are  specifically  designed  for  input-output  identification  experiments.  GSF  are 
based  on  information  theoretical  criteria  (the  Fisher  information  matrix)  and,  when  used  in 
conjunction  with  the  traditional  sensitivity  functions,  give  a  more  accurate  picture  on  the 
time  distribution  of  the  information  content  of  measured  outputs  with  respect  to  individual 
model  parameters. 

In  order  to  use  the  well-developed  sensitivity  analysis  for  the  theory  of  dynamical  systems 
for  our  purposes,  we  begin  by  writing  the  system  (9)  in  vector  form.  We  introduce  the 
notation 


c (t)  =  (c1(t),c2(t),c3(t),c4(t))T,  q  =  (ki,  ...,  K^Sm,l2,  -,U),  c0  =  (ci(0), ...,  c4(0))t, 


and  denote  by  F  =  (/i,  f2,  f3,  f±)T  the  vector  function  whose  entries  are  given  by  the  expres¬ 
sions  in  the  right  side  of  (9).  Then  F  :  l4  x  I8  ->  M4,  and  we  can  write  our  ODE  system  in 
the  general  vector  form 

dc 

Tt(t)  =  F(c,q),  (13) 

c(0)  =  c0. 

In  order  to  quantify  the  variation  in  the  state  variable  c (t)  with  respect  to  changes  in  the 
parameters  qj,  j  —  1, . . . , 8  and  the  initial  conditions  c0k,  k  =  1, ...  ,4,  we  are  naturally  led 
to  consider  the  sensitivity  matrices 


and 


3= b-,8 


Z  1=1,..., 4 

fc=l , ...  ,4 


dcj 

dqj 

dci 


dc{ 


Ok 


i=  1 , ...  ,4  ’ 

(14) 

1 ....  ,4 
k=l,.  ..’,4 

(15) 

We  note  that  since  our  function  F  is  sufficiently  regular,  the  solutions  ct  are  differentiable  with 
respect  to  qj  and  cok  and  therefore  our  sensitivity  matrices  Y  and  Z  are  well  defined.  The 
physical  interpretation  of  the  sensitivity  matrices  is  obvious.  Similar  to  the  partial  derivatives 
through  which  they  are  defined,  they  have  a  local  character  (in  time  and  parameters).  If, 
for  example,  the  entry  yi3  =  dci/dq3  of  the  matrix  Y  takes  values  very  close  to  zero  in  a 
certain  time  subinterval,  then  the  state  variable  c?;  is  insensitive  to  the  parameter  q3  on  that 
particular  subinterval.  The  same  entry  yi3  can  take  large  values  on  a  different  subinterval, 
indicating  that  in  this  time  subinterval,  the  state  variable  Cj  is  very  sensitive  to  the  parameter 
q3. 

From  sensitivity  analysis  theory  for  dynamical  systems  [10,  20,  25,  36,  42],  we  also  know 
that  Y (t)  is  a  4  x  8  matrix  that  satisfies  the  ODE  system 


Y(t)  =  Fc(c,  q)Y(t)  +  Fq(c,  q), 
Y(0)  —  04Xg, 


(16) 
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and  Z  (t)  is  a  4  x  4  matrix  that  satisfies 


Z (t)  =  Fc(c,  q)Z(t), 
Z(0)  =  I4x4. 


(17) 


Here  we  have  used  the  notation  Fc  =  <9F/<9c  and  Fq  =  dF/dq  for  the  4x4  and  the  4x8 
Jacobian  matrices  of  F  with  respect  to  c  and  q,  respectively,  while  0  and  I  are  the  zero  and 
the  identity  matrices  with  appropriate  dimensions.  Note  that  while  equations  (16),  (17)  are 
linear  in  Y  and  Z,  they  must  be  solved  in  tandem  with  equation  (13),  which  is  nonlinear. 
Consequently,  the  sensitivity  analysis  involves  the  solution  of  a  set  of  nonlinear  equations. 

We  will  compute  the  sensitivity  of  the  system  (13)  with  respect  to  q  and  Co  when  the 
solutions  are  essentially  at  steady  state.  We  carry  this  out  by  numerically  solving  the  systems 
(16)  and  (17)  for  the  same  values  of  the  parameters  q  =  q*  and  initial  conditions  c  =  cq  as 
used  in  the  stochastic  simulations  presented  above  (i.e. ,  those  given  in  Table  3)  and  by 
evaluating  the  solution  at  the  fixed  time  t  =  210  (arbitrarily  chosen,  but  sufficiently  large  for 
our  system  to  closely  approach  its  steady  state).  Due  to  the  nature  of  our  problem,  in  which 
the  parameters  have  different  units  and  the  state  variables  vary  widely  over  many  orders  of 
magnitude,  it  is  appropriate  to  consider  the  relative  sensitivities  SCiAj  defined  as  the  limit  of 
the  relative  change  in  c*  divided  by  the  relative  change  in  qj  when  the  relative  change  in  q 7- 
goes  to  zero,  i.e., 

A  a/ci 


S, 


■iiQj 


—  linr  .  , 

W  -o  A qj/qj 


(18) 


A  simple  analysis  of  the  definition  above  (assuming  that  both  c*  and  q3  are  nonzero)  yields 
that  the  relative  sensitivity  SCuqj  can  be  obtained  by  normalizing  the  usual  sensitivities 
dci/dqj  such  that 


o  _  m  <£ 


dcj 

dq, 


(19) 


We  note  that  the  SCi_q.  are  dimensionless  variables,  invariant  with  respect  to  changes  in  units 
for  Ci  and  qj,  which  we  can  utilize  to  compare  the  degree  of  sensitivity  of  the  state  variables 
with  respect  to  different  parameters.  In  Table  4,  we  tabulate  the  relative  sensitivities  at 
time  t  =  210  of  each  state  variable  c*  with  respect  to  each  parameter  qj  and  each  initial 
condition  Co*,.  For  any  fixed  parameter/initial  condition,  we  also  tabulate  the  sensitivity  of 
the  system ,  cumulatively  defined  as  the  Euclidean  norm  of  the  relative  sensitivities  of  the 
four  state  variables  with  respect  to  that  parameter /initial  condition.  In  other  words,  the 
sensitivity  of  the  system  with  respect  to  qj  is  given  by 


S„  = 


LA 


1/2 


(20) 


i=l 


For  the  particular  choice  of  the  parameters  q  =  q*  and  for  the  particular  initial  condition 
Co  =  Cq,  the  data  displayed  in  the  last  column  of  Table  4  reveals  that  near  the  steady-state 
(t  =  210)  the  system  (9)  is  most  sensitive  (in  decreasing  order)  to  l3,  c03  and  /2  and  least 
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Cl 

c2 

c3 

c4 

System 

K 1 

-2.65  x  10"1 

8.97  x  10~2 

8.34  x  10~3 

2.08  x  10“3 

2.80  x  10"1 

«  2 

-4.50  x  10”1 

-5.78  x  10”2 

8.76  x  10~2 

2.19  x  10~2 

4.62  x  10"1 

«3 

1.70  x  lO"1 

-9.05  x  10~3 

-2.40  x  10~2 

2.43  x  10"1 

2.98  x  10"1 

K4 

5.45  x  10"1 

-2.28  x  10"2 

-7.18  x  10-2 

-2.67  x  10”1 

6.12  x  10-1 

Sm 

0 

0 

0 

0 

0 

h 

-2.42  x  10° 

8.22  x  10"1 

7.65  x  10"2 

1.91  x  10"2 

2.56  x  10° 

h 

-5.17  x  10° 

-6.65  x  10"1 

1.00  x  10° 

2.51  x  10"1 

5.32  x  10° 

u 

6.80  x  10"1 

-3.61  x  10"2 

-9.63  x  10"2 

9.75  x  10"1 

1.19  x  10° 

Coi 

1.86  x  10° 

2.27  x  10"1 

2.11  x  10"2 

5.29  x  10"3 

1.88  x  10° 

C02 

1.63  x  10° 

1.99  x  10"1 

1.85  x  10"2 

4.63  x  10"3 

1.64  x  10° 

C03 

3.49  x  10° 

4.27  x  10"1 

3.97  x  10"2 

9.93  x  10"3 

3.52  x  10° 

0)4 

3.84  x  10”1 

4.69  x  10~2 

4.37  x  10~3 

1.09  x  10~3 

3.87  x  10"1 

Tabic  4:  Relative  sensitivities  of  the  state  variables/  system  with  respect  to  parameters  qj 
and  the  initial  conditions  c ok  near  steady  state  (t  =  210).  Baseline  parameter  values  (q  =  q*) 
and  initial  conditions  Cq  =  Cq  are  as  in  Table  3. 


sensitive  to  sm.  One  interesting  outcome  we  obtain  from  the  analysis  of  the  sensitivity  data 
presented  in  Tables  4  and  5  is  that  the  state  variable  ci  is  more  sensitive  to  the  parameter  I3 
and  the  initial  condition  C03  as  compared  to  the  parameters  ,  k4,  sm  and  /4  and  the  initial 
condition  C01  on  which  this  state  variable  depends  directly.  Without  this  sensitivity  analysis, 
we  could  not  infer  this  behavior  simply  by  looking  at  the  system  (9)  since  the  parameters 
I3  and  C03  do  not  appear  in  the  right  side  of  the  equation  which  defines  dc\jdt.  We  find 
similar  results  for  C2  and  c4  which  are  more  sensitive  to  the  non-direct  initial  condition  C03 
as  compared  to  their  direct  initial  conditions  C02  and  Co4. 

As  we  pointed  out  previously,  the  sensitivity  functions  by  definition  have  a  local  charac¬ 
ter  both  in  the  time  domain  and  in  the  parameter  domain,  which  implies  that  the  results 
displayed  in  Table  4  characterizes  the  sensitivity  of  the  system  only  for  q  =  q*,  Co  =  Cq  and 
t  =  210.  However,  to  obtain  a  broader  picture  of  the  sensitivity  map  for  our  system  (9)  near 
steady  state  ( t  =  210)  or  for  the  whole  time  interval,  one  can  compute  the  relative  sensitivi¬ 
ties  (19)  and  (20)  with  respect  to  q.j  and  Coa,  on  a  uniform  grid  in  a  parameter  neighborhood 
around  the  central  values  q*  and  c*Qk  and  analyze  the  results  thus  obtained.  (Although  we 
carried  out  such  analyses,  the  results  are  not  presented  here.) 

The  sensitivity  analysis  we  have  performed  so  far  is  usually  encountered  in  simulation 
studies  (direct  problems)  where  we  need  to  quantify  the  effects  of  parameter  variations 
011  the  trajectories  of  model  outputs.  Unlike  simulations,  in  identification  studies  (inverse 
problems)  one  typically  wants  to  estimate  model  parameters  from  data  measurements  and 
one  question  of  interest  is  to  determine  at  which  time  points  the  measurements  are  most 
informative  for  the  estimation  of  specific  parameters.  In  addition  one  may  also  desire  a 
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Cl 

c2 

c3 

C4 

System 

Sm 

Sm 

Sm 

Sm 

$171 

K-S 

«  3 

C04 

C04 

ftl 

ft  1 

ft4 

ftl 

ft  1 

^3 

C04 

k 

C02 

C02 

C04 

K  2 

C04 

C01 

C01 

«2 

ft4 

«3 

C03 

ft4 

k 

ft  1 

C03 

k 

k 

C02 

Co  2 

ft4 

«2 

C02 

Coi 

Coi 

k 

^3 

Coi 

k 

C03 

K2 

k 

k 

C03 

h 

k 

ft4 

C03 

h 

k 

h 

k 

k 

Table  5:  Summary  of  the  sensitivity  analysis  presented  in  Table  4.  The  columns  rank  the 
parameters  and  initial  conditions  in  order  of  increasing  sensitivity  for  the  four  state  variables, 
c*,  and  the  system. 


priori  information  about  the  degree  of  correlation  between  model  parameters  which  the 
TSF  functions  of  (14)  and  (15)  used  alone  cannot  provide.  To  address  these  questions, 
Thomaseth  and  Cobclli  [44]  introduced  the  generalized  sensitivity  functions  (GSF)  which 
provide  information  on  the  relevance  of  measurements  of  output  variables  of  a  system  for 
the  identification  of  certain  parameters  as  well  as  on  parameter  correlation.  More  precisely, 
the  generalized  sensitivity  functions  describe  the  sensitivity  of  the  parameter  estimates  with 
respect  to  data  measurements.  We  illustrate  below  the  utility  of  these  functions  to  provide 
a  better  understanding  of  our  network  model,  by  performing  a  sensitivity  analysis  for  the 
inverse  problem  of  estimating  the  parameters  ft’s  of  the  system  (9)  through  an  ordinary  least 
squares  procedure. 

For  a  single-output  model  f(t,9)  (e.g.,  the  case  where  one  has  longitudinal  observations 
of  one  component  c*  of  the  system  (13))  with  discrete  time  measurements 

y{h)  =  f(tk,90)  +ek,  k  =  1, . . .  ,Af,  (21) 


where  #o  is  the  “true”  parameter  values  (assumed  to  exist  in  most  statistical  inference  and 
information  content  formulations-see  [8,  11,  43]),  the  generalized  sensitivity  functions  (GSF) 
are  defined  by 


k  M  i 

gS  (4)  =  £{([£ 


i=l  j= 1 


-1  Vgf(ti,e0) 

X 


V„/(M 0)  ,  (22) 


where  ”  •”  represents  element-by-element  multiplication  (see  [44]  and  the  Appendix  of  [7]  for 
motivation  and  details).  The  measurement  errors  in  (21)  are  assumed  to  be  independently 
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and  identically  distributed  with  zero  mean  and  known  variance  a 2(fy.)  and  the  nonlinear 
model  function  fit,  6)  is  assumed  to  be  differentiable  with  respect  to  9.  Moreover,  for 
simplicity,  it  is  assumed  that  the  observations  y(tk),  as  represented  in  (21),  are  used  to 
estimate  90  by  minimizing  the  weighted  sum  of  squares  (WSS) 


WSS(y,9 ) 


V  fa/fc)  -  /(M)]2 

rr 


(23) 


In  our  case,  the  nonlinear  model  function  /  is  replaced  by  a  vector-valued  function  which  is 
the  solution  of  the  system  (9)  and  the  generalized  sensitivity  functions  (GSF)  are  given  by 


k  4 

gs(tfc)  =  EE{([EE  a2u  yeCi(tjl90)Vdci(tjl90)': 


M 


i=  1  1=1 


= l  i=i  1  y 


1  VgCi(ti,90) 
X 


tfiti) 


VgCi(ti,90) 


(24) 

Here  VgCi  represents  the  gradient  of  the  state  variable  q  with  respect  to  9 ,  where  9  €  Mp  is  a 
vector  including  all  (or  just  a  subset  of)  the  parameters  Ac’s  and  Z’s  and  the  initial  conditions 


c0’s. 


We  note  that  the  generalized  sensitivity  functions  (24)  are  vector-valued  functions  having 
the  same  dimension  P  as  9  and  defined  only  at  the  discrete  time  points  tk,k  =  1, . . . ,  M.  They 
are  cumulative  functions,  at  each  time  point  tg  taking  into  account  only  the  contributions  of 
the  measurements  up  to  A,  thus  representing  the  influence  of  longitudinal  measurements  in 
contributing  to  the  parameter  estimates. 

From  (24)  it  follows  that  all  the  components  of  gs  are  one  at  the  end  of  the  experiment, 
i.e. ,  gs (£m)  =  1.  If  one  defines  gs (t)  —  0  for  t  <  ti  (gs  is  zero  when  no  measurement  is 
collected)  and  interpolates  gs  continuously  between  observation  times,  then  each  component 
gsp  of  gs  varies  continuously  from  0  to  1  during  the  experiment.  As  we  will  see  in  the  example 
below,  this  transition  is  not  necessarily  monotonic  ( gsp,p  —  1, . . . ,  P  may  have  oscillations) 
nor  is  it  restricted  to  values  in  [0, 1]  (i.e.,  gsp  may  take  values  outside  [0, 1])  if  large  corre¬ 
lations  between  parameter  estimates  exist.  As  discussed  in  [44],  the  time  subinterval  during 
which  this  transition  has  the  sharpest  increase  corresponds  to  measurements  which  provide 
the  most  information  on  possible  variations  in  the  corresponding  true  model  parameters. 

Since  the  GSF  theory  is  developed  in  the  context  of  estimation  problems,  we  considered 
next  an  estimation  problem  using  simulated  “data”.  The  data  used  for  inversion  was  gen¬ 
erated  by  first  numerically  solving  the  system  (9)  for  the  parameter  values  given  in  Table  3 
and  then  adding  5%  Gaussian  white  noise  to  the  solution  obtained.  We  consider  the  problem 
of  using  this  data  to  estimate  the  parameters  ACi,  AC2,  AC3  and  AC4  in  an  ordinary  least  squares 
procedure  (with  the  other  parameters  and  initial  conditions  fixed  at  the  values  from  Table 
3).  The  true  values  for  the  Ac’s  are  R  =  (1.674,0.322,4.521, 1)T.  For  90  =  R,  the  generalized 
sensitivity  functions  (24)  along  with  the  traditional  sensitivity  functions  for  the  system  (9) 
are  presented  in  Figure  7.  I11  both  figures  we  note  a  very  well  defined  time  subinterval,  from 
t  —  0  to  about  t  =  60,  where  both  GSF  and  TSF  plots  exhibit  sharp  increases/decreases. 
After  this,  the  TSFs  reach  very  quickly  a  steady  state  and  the  GSFs  are  forced  to  approach 
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Generalized  Sensitivity  Functions  with  respect  to  k2,k3,k4 


Figure  7:  Generalized  and  Traditional  Sensitivity  Functions  for  K\,  k2,  k3,  k4,  employing  a 
simulated  data  set  with  M  =  210,  generated  as  discussed  in  the  text.  Underlying  parameter 
values  are  K\  =  1.674,  k2  =  0.322,  k4  =  4.521,  k4  =  1,  while  all  other  parameters  are  as 
given  in  Table  3. 


one.  According  to  the  theory,  the  interval  [0,60]  is  the  time  region  where  measurements 
are  most  informative  for  estimating  the  true  parameters  R.  So  at  least  intuitively,  sampling 
more  data  points  in  this  region  would  result  in  more  information  about  the  parameters  R 
and  therefore  more  accurate  estimates  for  them.  By  computing  the  correlation  matrix  whose 
elements  are  given  by  standard  formulas  in  least  squares  theory  [9,  43],  one  can  also  observe 
that  strong  correlations  exist  between  estimates  for  k3  and  k4.  In  fact,  the  correlation  matrix 
for  these  parameters  is  given  by 


Corr 

«i 

re  3 

k4 

K 1 

1.0000 

0.4941 

0.0004 

0.2209 

re  2 

0.4941 

1.0000 

0.1388 

0.1497 

0.0004 

0.1388 

1.0000 

-0.9502 

K4 

0.2209 

0.1497 

-0.9502 

1.0000 

which  is  in  agreement  with  the  dynamics  of  the  curves  shown  in  Figure  7.  Positive  correlation 
between  K\  and  k2  is  clearly  indicated  as  the  corresponding  gsf  graphs  increase  together  while 
the  negative  correlation  between  k3  and  k4  is  evidenced  by  the  early  opposite  slope  behavior 
in  their  corresponding  gsf  graphs. 

Ordinary  least  squares  inverse  problems  carried  out  with  different  sets  of  data  points 
illustrate  and  support  the  theory.  We  first  performed  the  least  squares  minimization  for  a 
data  set  DSi  consisting  of  a  total  of  15  observations,  of  which  8  were  taken  at  equidistant 
points  in  the  interval  [0,  60]  and  7  taken  at  equidistant  points  in  the  interval  [80,  210]  (for 
simplicity,  we  exclude  the  transition  interval  [60,80]  from  our  analysis  here).  The  estimates 
for  Ki,  k2,  k3  and  k4  along  with  the  M4  Euclidian  norm  of  the  error  are  displayed  in  Table  6. 
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Data  Points 

in 

Estimates  for 

Data  Set 

[0,60] 

[80,210] 

Total 

K.\ 

K  2 

K4 

error 

DSi 

8 

7 

15 

1.638 

0.329 

1.775 

2.028 

2.932 

ds2 

15 

7 

22 

1.579 

0.305 

4.341 

0.966 

0.207 

ds3 

30 

7 

37 

1.708 

0.302 

2.728 

1.178 

1.803 

ds4 

8 

14 

22 

1.625 

0.322 

1.830 

1.786 

2.804 

ds5 

8 

0 

8 

1.686 

0.337 

2.379 

1.597 

2.224 

ds6 

15 

0 

15 

1.736 

0.330 

6.378 

1.064 

1.859 

ds7 

30 

0 

30 

1.724 

0.298 

3.549 

1.013 

0.973 

DSs 

0 

17 

17 

12.540 

1.966 

14.587 

9.814 

17.314 

DSg 

0 

33 

33 

45.554 

7.588 

69.183 

32.589 

84.600 

DSW 

0 

66 

66 

51.161 

9.842 

92.137 

36.002 

106.964 

True  Parameter  Values 

1.674 

0.322 

4.521 

1 

Table  6:  Parameter  estimates  for  the  transition  rates  K\,  k2 ,  k3  and  K4  with  different  data 
sets. 


If  we  increase  the  number  of  data  points  in  the  interval  [0,  60]  from  8  to  15  and  keep  the 
number  of  data  points  in  [80,210]  the  same  (see  DS2  entry,  same  table),  we  observe  a 
significant  decrease  in  the  Euclidian  norm  of  the  error  from  2.932  to  0.207  which  represents 
a  significant  increase  in  the  accuracy  of  the  parameter  estimates.  A  similar  decrease  from 
2.224  to  1.859  and  then  to  0.973  is  observed  when  we  solve  the  least  squares  problem  only 
with  data  from  [0,  60]  (see  DS$,  DSq  and  DS7  entries).  Thus,  numerical  calculations  support 
the  fact  that  increasing  the  number  of  data  points  in  the  region  [0,60]  yields  more  accurate 
estimates  for  the  parameter  R,  in  agreement  with  the  theoretical  expectations  from  TSF  and 
GSF. 

A  totally  different  outcome  is  obtained  when  we  carry  out  numerical  estimation  with  an 
increasing  number  of  data  points  in  the  interval  [80,210].  As  one  can  see  by  comparing  the 
entries  DS\  and  DS4  in  Table  6,  only  a  small  gain  is  obtained  in  the  accuracy  of  the  parameter 
estimates  is  gained  by  doubling  the  number  of  data  points  in  the  interval  [80,  210].  Moreover, 
when  we  try  to  estimate  the  parameters  R  with  data  from  the  interval  [80,  210]  alone  (see 
DSg ,  DSg  and  DS\0),  we  obtain  very  large  errors  which  increase  in  magnitude  as  the  number 
of  sample  points  increases.  Although  puzzling  at  first  view,  this  phenomenon  is  not  surprising 
at  all  from  the  perspective  of  the  theory  presented  above  and  that  presented  in  [6] .  Indeed,  by 
blindly  sampling  more  data  points  from  the  region  where  the  generalized  sensitivity  functions 
exhibit  the  so  called  “forced-to-one”  artifact  and  the  traditional  sensitivity  curves  are  flat,  we 
simply  introduce  redundancy  in  the  sensitivity  matrix,  thus  increasing  the  condition  number 
of  the  Fisher  information  matrix  for  our  problem.  For  an  illustration  and  discussion  of  this 
phenomenon,  see  [6].  By  the  Cramer- Rao  inequality,  the  consequence  is  that  the  variance 
of  the  unbiased  estimator  (and  the  corresponding  standard  errors)  will  be  huge,  making  our 
estimates  less  useful.  The  same  phenomenon  (introduction  of  redundancy  in  the  sensitivity 
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matrix)  is  responsible  for  the  poorer  results  which  are  obtained  when  we  estimate  R  using 
data  set  DS%  (the  Euclidian  norm  of  the  error  increased  as  compared  to  DS2,  where  we 
doubled  the  number  of  points  in  [0,60]).  We  observe  an  important  drawback  of  the  GSF: 
while  they  specify  the  most  informative  regions  with  respect  to  the  estimation  of  parameters, 
they  do  not  provide  any  information  about  the  necessary  number  of  data  points  to  be  used 
in  those  regions. 

We  remark  that  in  this  section  we  have  illustrated  a  methodology  to  quantify  the  sen¬ 
sitivity  with  respect  to  parameters  and  initial  conditions  of  a  large  complex  system.  This 
methodology  is  a  fundamental  tool  in  identifying  those  parameters  in  a  model  which  require 
further,  more  in-depth  investigation  as  well  as  the  basis  of  a  statistical  analysis  for  quantify¬ 
ing  uncertainty  in  estimators  (e.g.,  see  [5,  6,  9,  11,  43])  as  well  as  information  content  based 
model  selection  techniques  [8]. 

4  Foot-and-Mouth  Disease 

Having  established  a  basic  model  for  the  movement  of  animals  (pigs)  in  the  agricultural 
system  (the  pork  industry  of  North  Carolina),  we  are  now  ready  to  model  the  spread  of  an 
infectious  disease  in  the  food  production  network.  In  this  section  we  describe  the  incorpora¬ 
tion  of  an  SIR-type  infection  into  the  system  and  present  simulations  to  illustrate  the  spread 
of  Foot  and  Mouth  Disease  throughout  the  aggregated  agricultural  network. 

We  describe  the  infection  by  an  SIR  process  [2,  12].  It  is  assumed  that  a  population  can 
be  partitioned  into  three  groups:  susceptible  (S),  infectious  (I)  and  removed  (R).  In  many 
settings  the  removed  class  represents  individuals  who  have  recovered  from  the  infection. 
Individuals  move  between  these  classes  as  they  become  infected  and  recover  from  infection. 
Recovery  is  assumed  to  confer  permanent  immunity  to  infection  and  the  demography  of  the 
population  (i.e. ,  births  and  non  disease-related  deaths)  is  ignored.  In  a  well-mixed  population 
the  epidemiological  model  can  be  described  by  the  flowchart  of  Figure  8  and  equations  (25). 

Susceptible  Infectious  Recovered 


Figure  8:  Flow  diagram  of  the  simple  SIR  model. 
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Here,  S,  I  and  R  denote  the  numbers  of  susceptible,  infective  and  removed  individuals, 
respectively.  The  transmission  parameter  is  (3:  this  parameter,  when  combined  with  the  rate 
at  which  individuals  meet  each  other  and  the  probability  that  an  infective  would  infect  a 
susceptible  during  any  one  such  meeting,  yields  the  transmission  rate.  It  is  assumed  that 
recovery  occurs  at  constant  rate  7,  so  that  I/7  is  the  average  duration  of  infection.  The 
population  size  is  denoted  by  N,  and  we  have  in  this  case  that  N  —  S  +  I  +  R. 

The  behavior  of  the  simple  SIR  model  is  governed  by  the  basic  reproductive  number, 
R0  [2,  12,  31,  14].  This  quantity  equals  the  number  of  secondary  infections  caused  by 
the  introduction  of  a  single  infectious  individual  into  an  otherwise  completely  susceptible 
population.  In  terms  of  model  parameters,  the  basic  reproductive  number  is  given  by 

Ro  =  P/r 

An  outbreak  of  infection  can  only  ever  occur  if  Ro  is  greater  than  one,  otherwise  the  number 
of  infectives  can  never  increase. 

We  now  combine  the  agricultural  network  model  and  the  SIR  model  to  produce  a  de¬ 
scription  of  the  potential  spread  of  an  SIR-type  infection  in  our  agricultural  system.  It  is 
most  convenient  for  us  to  work  with  the  numbers  of  individuals  of  each  type  found  at  the 
various  nodes  of  the  network,  and  so  we  convert  the  concentrations  of  animals  of  the  network 
model  (9)  into  numbers.  We  write  the  number  of  individuals  found  at  node  i  as  Ni}  and  so 
we  have  that  Ni(t)  =  iVcj(f). 

We  expand  the  first  three  nodes  of  the  network,  i.e.,  those  describing  the  growers/sows, 
nurseries  and  finishers,  by  including  an  SIR  model  within  each  of  them.  The  infection 
statuses  of  the  animals  at  the  All  node  are  tracked  by  the  quantities  Si,  It  and  Ri.  Since  N,t 
is  defined  to  be  the  total  number  of  animals  at  this  node  we  have  that  N,  —  St  +  A  +  R,t. 

We  make  the  following  set  of  assumptions  about  the  infection  process  and  its  impact 
upon  the  agricultural  system: 

(1)  Pigs  are  born  healthy,  but  susceptible.  Piglets  are  introduced  into  the  network  at  the 
first  node. 

(2)  There  is  no  infection  or  recovery  during  transport  between  the  nodes.  This  is  based 
on  the  idea  that,  in  most  cases,  transportation  takes  no  more  than  a  couple  of  days, 
which  is  relatively  short  compared  to  the  amount  of  time  the  pigs  spend  at  each  node. 

(3)  There  is  no  infection  in  the  slaughter  node.  I11  the  absence  of  human  intervention  this 
does  not  mean  that  the  infectious  animals  are  not  processed:  the  assumption  is  that 
the  infection  is  not  propagated  after  the  animal’s  death. 

(4)  Infected  pigs  in  node  i  recover  at  rate  7*.  These  recovery  rates  may  differ  between  the 
nodes  and,  since  the  individuals  found  at  different  nodes  will  have  different  ages,  these 
parameters  depict  age-dependent  recovery  rates. 

(5)  Recovered  pigs  have  temporary  immunity,  i.e.,  pigs  do  not  immediately  become  suscep¬ 
tible  after  recovering  from  FMD.  The  rates  at  which  recovered  pigs  become  susceptible 
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again  are  pi  (i  —  1,  2  or  3).  While  some  diseases  afford  permanent  immunity  to  the  re¬ 
covered  animal,  immunity  acquired  after  FMD  infection  wanes  in  a  matter  of  months. 
If  vaccination  of  animals  were  considered  there  would  be  an  accumulation  of  animals 
in  the  recovered  groups  of  the  appropriate  nodes. 

(6)  We  assume  that  the  FMD-related  death  rate  is  small  enough  that  we  can  ignore  such 
deaths.  Control  strategies  such  as  culling  infective  or  susceptible  animals  could  be 
modeled  by  including  death  terms. 

(7)  Since  our  network  model  assumes  that  the  system  is  closed,  we  assume  that  any  deaths 
are  replenished  by  the  introduction  of  piglets  into  the  network.  As  already  men¬ 
tioned,  these  introductions  occur  into  the  first  node.  Consequently,  as  illustrated  in 
the  flowchart,  Figure  9,  the  slaughtered  animals  that  leave  node  four  effectively  return 
to  the  Si  class.  Similarly,  if  animal  culling  were  being  considered,  the  model  would 
include  flows  from  the  appropriate  classes  back  to  Si. 

(8)  No  human  intervention.  In  the  model  as  presented  here,  we  assume  that  humans  do 
not  make  any  adjustments  to  operation  of  the  agricultural  system  in  response  to  the 
infection:  animal  movement  and  processing  continues  as  normal.  Of  course,  one  of  the 
main  aims  of  creating  a  model  such  as  this  is  to  enable  the  consideration  of  control 
strategies.  In  this  study,  we  do  not  do  this  as  we  wish  to  first  establish  the  baseline 
(no-control)  behavior  of  the  system. 

Grower  Nursery  Finisher  Slaughter 


Figure  9:  Flow  diagram  of  the  aggregated  agricultural  network  model  with  SIR  infection. 

With  this  set  of  assumptions  we  obtain  the  following  system  of  equations  for  the  deter¬ 
ministic  model  of  the  aggregated  agricultural  network  model  with  an  SIR  infection. 


26 


7  Q 

~ jj~  —  —fliSiIi/Ni  +  p1R1  —  kiSi(L2  —  N2)+  +  h  min(A^4,  Sm ax) 

~ jj ■  =  fliSih/Ni  —  '-fill  —  k\I\{L2  —  N2)+ 

dU 

— =  7iii  —  pi-Ri  —  k\R\{L2  —  N2)+ 

7  Qf 

=  —(32S2I2/N2 -\-  p2R2  —  k2S2(Lz  —  Nz)+ -\- k\S\(L2  —  N2)+ 

—jj-  =  fd2  S2 1 2  /  N2  —  y2I2  —  k2I2(L3  —  N3)+  +  k\I\{L2  —  N2)+ 

dR 

-jjj-  =  72-^2  —  P2R2  ~  k2R2{L3  —  N3)+  +  k\R\(L2  —  N2)+ 

7  Qf 

-jjj-  =  —  (33S3I3/N3  +  P3R3  —  k3S3(L4  —  N4)+  +  k2S2(L3  —  ^3)4. 

dTn 

-jj-  =  P3S3I3/N3  —  73/3  —  k3I3(L4  —  N4)+  +  k2I2(L3  —  N3)+ 

dR- 

-jjj-  =  73-^3  —  P3-R3  —  k3R3{L4  —  N4)+  +  k2R2{L3  —  N3)+ 

—jjj—  =  min(7V4,  5max)  +  k3N3(L4  -  N4)+. 

In  many  ways,  this  model  resembles  a  standard  multi-group  epidemiological  model  and 
so  we  might  hope  to  be  able  to  find  the  basic  reproductive  number  of  the  system  using 
standard  multi-group  methodology  [12,  13].  It  is  straightforward  to  find  the  next  generation 
matrix,  whose  entries,  tij,  give  the  average  numbers  of  secondary  infections  that  result  in 
node  i  from  the  introduction  of  one  infectious  individual  into  node  j  when  all  individuals  at 
node  %  are  susceptible.  For  within-patch  transmission,  we  have 

_ As _ 

73  +  k3(L4  —  IV4) 

_ Z-^2 _ 

72  +  k2{L3  —  N3)+ 

Pi 

7i  +  k\(L2  —  N2)+ 

Here,  the  term  1/(7*  +  ki(Li+ 1  —  fV*+i)+)  is  the  average  duration  of  infection  for  infectives  in 
node  i,  corrected  for  their  departure  on  account  of  transportation.  Movement  of  individuals 
between  nodes  reduces  the  average  number  of  within-node  secondary  infections,  with  this 
effect  being  most  noticeable  if  the  node  residence  time  is  short  compared  to  the  duration  of 


£33  = 
^22  = 
tn  = 
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infection.  For  between-node  transmission  we  have 

^2(^3  —  Af3)+  , 

72  +  k2(L3  —  ^3)+ 
k\{L2  —  N2)+ 

7i  +  ki(L2  —  N2)+  22 

ki(L2  —  N2)+  k2(L3  —  N3)+  ^ 

7i  +  ki(L2  —  N2)+  72  +  k2(L3  —  N3)+ 

^23  —  ^13  —  0. 

Here,  the  term  ki(Li+ 1  —  7Vj+i)+/(7*  +  ki(Li+ 1  —  7Vi+i)+)  gives  the  probability  that  an  infec¬ 
tious  individual  in  node  i  is  transported  to  node  i  +  1  before  it  recovers. 

Multi-group  theory  reveals  that,  for  an  irreducible  system  (i.e.,  one  in  which  infection  can 
travel  between  any  pair  of  nodes,  possibly  via  intermediate  nodes),  the  basic  reproductive 
number  is  given  by  the  dominant  eigenvalue  of  this  matrix.  Our  system  is  not  irreducible, 
however,  since  infected  animals  can  only  move  from  a  node  to  the  following  node:  infection 
can  spread  to  succeeding  nodes  but  not  preceding  ones.  Consequently,  the  standard  definition 
of  Rq  for  multi-group  systems  is  not  helpful  to  us.  From  the  tij  we  can,  however,  easily 
calculate  the  average  number  of  secondary  infections  caused  by  the  introduction  of  a  single 
infective  into  node  j  when  all  other  individuals  in  the  population  are  susceptible.  These 
quantities  equal  tn  +  t23  +  t3 1,  t22  +  t32  and  t33,  for  nodes  1,  2  and  3,  respectively. 

We  now  turn  to  numerical  simulation  of  the  model.  Unfortunately,  the  existing  stud¬ 
ies  in  the  literature  do  not  provide  us  with  an  appropriate  set  of  parameter  values  to  use. 
Epidemiological  parameters  for  the  spread  of  FMD  between  various  types  of  animals,  in¬ 
cluding  pigs,  have  been  quantified  [17,  18,  22,  23,  24,  32,  33,  34,  38,  39,  46]  but  on  spatial 
scales  that  are  quite  different  from  our  aggregated  network  description.  On  a  large  spatial 
scale,  transmission  between  farms  has  been  described,  for  instance  during  the  1967/68  and 
2001  outbreaks  in  the  UK,  and  parameter  estimates  obtained  [22,  23,  24,  32,  33,  34,  46]. 
Large-scale  studies,  however,  take  the  individual  unit  of  the  model  to  be  farms  and  so  do 
not  consider  transmission  between  individual  animals. 

On  a  small  spatial  scale,  detailed  transmission  experiments  [17,  18,  39,  38]  have  examined 
the  spread  of  infection  between  small  numbers  of  closely  housed  animals,  either  within  or 
between  pens.  (Given  the  earlier  comment  regarding  age-dependent  epidemiological  param¬ 
eters,  it  is  interesting  to  note  that  some  of  these  experiments  have  been  carried  out  on  pigs 
of  different  ages.)  These  experiments,  which  typically  involve  placing  one  or  more  infected 
animals  in  close  contact  with  a  number  of  susceptible  animals,  demonstrate  the  high  degree 
of  infectiousness  of  FMD.  In  many  instances  every  susceptible  animal  became  infected  [17]. 
Instances  in  which  all  animals  become  infected  provide  less  informative  estimates  of  R0  than 
might  be  hoped  since  the  statistical  methodology  employed  gives  an  infinite  estimate  for  R0. 
An  alternative  statistical  approach  [18]  accounts  for  the  time  dependence  in  the  experiment 
and  provides  estimates  of  the  transmission  parameter  (3 . 

The  literature  provides  us  with  satisfactory  estimates  for  the  average  durations  of  infec¬ 
tiousness  and  immunity  [28],  but  not  R0  (or,  equivalently,  (3).  For  illustrative  purposes,  we 


£32  = 

£21  = 

£31  = 
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shall  take  Rq  to  equal  10  at  each  level  of  the  network  in  the  absence  of  transportation. 


Parameters  for  the  Aggregated  Agricultural  Network  Model  with  SIR 

Parameters 

Definition 

Values 

Units 

Reference 

Rq 

basic  reproductive  number 

10 

7u72,73 

recovery  rate 

1/31,1/31,1/31 

1/  day 

[28] 

PliP2,P3 

rate  of  loss  of  immunity 

1/180,1/180,1/180 

1/  day 

[28] 

All  the  simulations  that  follow  assume  that  the  network  is  initially  in  its  demographic 
steady  state.  The  epidemiological  assumptions  that  we  have  made  lead  to  this  holding  for 
all  future  times.  (This  would  not  be  the  case  if  we  had  deaths  at  the  nursery  or  finisher 
nodes.)  We  choose  to  introduce  infection  by  infecting  a  certain  number  of  individuals  in  the 
initial  node. 

Our  first  simulation  (Figure  10)  illustrates  the  movement  of  a  cohort  of  infectious  individ¬ 
uals  through  the  network  in  the  absence  of  ongoing  transmission.  We  set  /A  =  /?2  =  @3  =  0, 
and  introduce  infection  by  placing  50,  000  infected  piglets  in  the  grower  node. 

In  this  simulation  infected  animals  either  recover  or  get  transported  to  the  next  node. 
In  the  grower  node,  the  number  of  infected  animals  simply  decreases  exponentially,  while 
the  number  of  recovered  animals  increases  and  then  declines.  Within  a  fairly  short  time 
period,  the  grower  population  is  entirely  replaced  by  susceptible  individuals,  reflecting  the 
rapid  turnover  of  the  grower  node.  We  see  the  appearance  of  infection  first  at  the  nursery 
and  then  at  the  finisher  node.  Since  the  system  is  at  demographic  equilibrium,  we  see  no 
change  at  the  slaughter  node. 

For  our  second  simulation  (Figure  11),  we  assume  that,  in  the  absence  of  transportation, 
the  infection  would  have  basic  reproductive  number  equal  to  10  at  each  node.  Consequently, 
we  set  Pi/'yt  =  10.  As  before,  we  take  the  initial  population  to  be  in  demographic  equilibrium, 
but  we  now  introduce  just  200  infective  piglets  into  the  grower  node. 

The  system  rapidly  approaches  an  equilibrium  in  which  infectious  animals  are  present  at 
each  node.  The  fraction  of  animals  that  are  susceptible  at  equilibrium  is  much  higher  at  the 
grower  node  (24.8%)  than  at  either  the  nursery  (6.9%)  or  finisher  (7.9%).  This  should  be 
expected  since  all  animals  arriving  at  the  grower  are  susceptible,  while  the  infection  statuses 
of  those  entering  the  nursery  or  finisher  reflect  the  composition  of  the  preceding  node  (i.e., 
grower  and  nursery  nodes,  respectively).  For  example,  only  25%  of  the  arrivees  at  the  nursery 
are  susceptible.  The  susceptible  fraction  at  the  finisher  is  slightly  higher  than  at  the  nursery 
due  to  the  replenishment  of  susceptibles  by  loss  of  immunity:  this  relatively  slow  process  has 
a  higher  chance  of  occurring  at  the  finisher  since  an  animal’s  average  residence  time  there  is 
longer  than  at  the  nursery. 

The  differing  susceptible  fractions  between  the  nodes  mean  that,  relative  to  the  size  of  the 
population  of  each  node,  there  is  less  ongoing  transmission  at  either  the  nursery  or  finisher 
node  than  at  the  grower.  The  equilibrium  infective  fraction  decreases  as  the  supply  chain  is 
traversed  (46.8%,  31.5%  and  16.2%  at  the  grower,  nursery  and  finisher,  respectively),  while 
the  recovered  fraction  increases  (28.4%,  61.5%  and  75.9%). 
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Figure  10:  Simulation  I:  Recovery  and  movement  of  infectious  individuals  through  the  net¬ 
work,  in  the  absence  of  transmission.  Notice  that  susceptible  numbers  are  plotted  on  different 
scales  (left  vertical  axis  on  each  panel)  from  infective  and  recovered  numbers  (right  vertical 
axis  on  each  panel).  Parameter  values  are  given  in  the  text:  each  (3i  is  set  equal  to  zero, 
infection  is  assumed  to  last  31  days  and  immunity  lasts  180  days.  The  initial  conditions  are 
Ni  =  315,  000,  N2  =  735,  000,  N3  =  2, 100,  000  and  N4  =  15, 000.  Of  the  piglets  at  the  first 
node,  50,  000  are  initially  taken  to  be  infectious,  while  the  remainder  of  the  population  is 
susceptible. 


The  impact  of  age-dependent  transmission  rates  and  the  differing  residence  times  is  ex¬ 
plored  in  our  third  simulation  (Figure  12).  Again  only  200  infectives  were  introduced  but 
now  we  assume  that  younger  animals  are  less  infectious  than  older  animals.  The  basic  repro¬ 
ductive  numbers,  in  the  absence  of  transportation,  at  the  grower,  nursery  and  finisher  nodes 
are  taken  to  equal  two,  ten  and  fifteen,  respectively.  (We  exaggerate  the  age-dependent 
differences  in  R0  for  illustrative  purposes.) 

Transportation  has  a  major  impact  on  the  dynamics  of  the  infection:  the  movement  of 
individuals  out  of  the  grower  node  is  sufficiently  rapid  that  the  infection  cannot  persist  in  this 
node.  The  number  of  infectives  falls  roughly  exponentially,  as  does  the  number  of  recovereds, 
although  the  latter  only  occurs  after  an  initial  rise  (which  largely  reflects  the  recovery  of  the 
initial  pool  of  infectives).  The  dynamics  in  this  node  are  somewhat  reminiscent  of  those  seen 
in  the  first  simulation  (in  which  there  was  no  transmission),  although  they  play  out  on  a 
slower  timescale  since  there  is  some  ongoing  transmission. 
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Figure  11:  Simulation  II:  Dynamics  following  the  introduction  of  infection,  with  i?o  =  10. 
This  value  of  the  basic  reproductive  number  determines  the  value  of  the  transmission  pa¬ 
rameters  (/3i).  Initial  conditions  and  other  parameter  values  are  as  in  the  previous  figure, 
except  that  only  200  of  the  315,000  individuals  at  the  grower  node  are  infective  at  the  initial 
time. 

Disease  transmission  at  the  nursery  and  finisher  nodes  occurs  sufficiently  rapidly  that 
the  prevalence  of  infection  approaches  a  positive,  endemic,  equilibrium  at  both.  Observe 
that,  even  though  the  transmission  parameter  in  the  nursery  node  is  the  same  as  it  was  in 
the  previous  simulation,  the  equilibrium  numbers  of  susceptibles  and  infectives  are  higher 
here  than  they  were  in  the  previous  simulation.  This  reflects  the  differences  between  the 
compositions  of  the  populations  entering  the  nursery  node  in  the  two  simulations,  with  more 
susceptibles  arriving  in  the  age-dependent  situation. 
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Figure  12:  Simulation  III.  Age-dependent  transmission  of  infection.  If  there  were  no  trans¬ 
portation  of  animals,  the  basic  reproductive  numbers  at  the  grower,  nursery  and  finisher 
nodes  would  equal  2,  10  and  15,  respectively.  These  values  of  Ro  determine  the  parameters 
j3i,  fa  and  /?3.  All  other  parameter  values  and  initial  conditions  are  as  in  the  previous  figure. 
For  the  panel  depicting  the  grower,  note  that  susceptible  numbers  are  plotted  on  a  different 
scale  (left  vertical  axis)  from  the  infective  and  recovered  numbers  (right  vertical  axis),  and 
these  infective  and  recovered  numbers  are  not  in  thousands.  Also  note  that  infection  goes 
extinct  in  the  grower  node  while  a  positive  equilibrium  level  of  infection  is  achieved  at  the 
nursery  and  finisher  nodes. 
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5  Concluding  Remarks 

In  this  paper  we  have  demonstrated  a  methodological  approach  to  the  investigation  of  pro¬ 
duction  networks  and  their  vulnerability  to  disturbances  such  as  diseases.  The  stochastic 
model  and  the  resulting  approximate  deterministic  system  we  employ  were  shown  to  agree 
well,  but  are  not  validated.  Rather,  we  carry  out  simulations  and  sensitivity  analyses  with 
parameter  values  that  are  only  loosely  based  on  a  swine  network.  We  use  the  deterministic 
model  to  show  how  to  determine  the  parameters  to  which  the  model  at  these  parameter 
values  exhibits  the  most  sensitivity.  Finally,  we  demonstrate  how  disease  can  be  introduced 
and  the  resulting  network  vulnerabilities  analyzed.  An  interesting  next  step  would  involve 
obtaining  experimental  data  to  validate  and  perhaps  improve  the  model  for  a  specific  pro¬ 
duction  network.  This  would  require  using  inverse  problem  algorithms  with  the  data  to 
obtain  estimates  along  with  measures  of  associated  uncertainty  (e.g.,  standard  errors  [9,  43]) 
for  the  underlying  transition  rates  kt. 

Other  obvious  questions  for  further  investigation  involve  the  introduction  of  stochasticity 
in  the  transmission  dynamics  of  the  infection.  It  is  well  known,  for  example,  that  random 
effects  can  have  a  major  impact  on  the  invasion  of  an  infection  into  a  population.  The  use  of 
constant  rates  of  recovery  and  loss  of  immunity  for  the  infection  should  also  be  questioned. 
These  assumptions,  which  are  biologically  unrealistic  for  many  infections,  could  be  important 
in  cases  such  as  the  one  presented  above  where  the  life  span  of  the  animals  is  comparable  in 
length  to  the  duration  of  infection  and  immunity. 

The  randomness  seen  in  the  stochastic  network  model  originates  from  the  random  move¬ 
ment  of  discrete  individuals  from  node  to  node.  The  analysis  of  Sections  2.2  and  3.2  shows 
that,  due  to  an  averaging  effect,  these  random  effects  become  less  important  as  the  system 
size  N  increases.  Application  of  the  stochastic  transportation  model  to  describe  a  real-world 
situation  should,  therefore,  account  for  the  size  of  the  groups  in  which  pigs  are  transported 
between  nodes.  If,  for  example,  one  thousand  pigs  were  moved  at  a  time,  the  appropriate 
notion  of  an  “individual”  within  the  model  would  be  a  thousand  pigs.  Treating  each  group 
of  a  thousand  animals  as  a  unit  would  lead  to  a  marked  increase  in  the  magnitude  of  stochas¬ 
tic  fluctuations  seen  at  the  population  level.  Consequently,  even  though  the  system  size  in 
our  simulations  is  on  the  order  of  millions  of  pigs,  it  might  be  that  the  resulting  stochastic 
fluctuations  in  a  more  realistic  model  of  the  production  system  are  closer  to  those  shown  in 
Figure  3  than  to  those  of  Figure  2. 

The  approach  outlined  in  this  paper  has  rather  obvious  potential  for  application  to  a 
wide  range  of  problems.  These  include  the  investigation  of  the  spread  of  diseases  through 
spatially  or  structurally  distributed  dynamic  populations  (e.g.,  avian  flu  through  migrating 
bird  populations,  contagious  infections  through  highly  mobile  and/or  age-structured  human 
and  animal  populations).  In  some  of  these  cases  the  natural  nodal  structure  would  be  a  con¬ 
tinuum,  requiring  stochastic  and  deterministic  models  with  a  continuum  of  spatial/structural 
heterogeneities,  leading  to  partial  differential  equation  systems.  Such  applications  would  un¬ 
doubtedly  motivate  the  development  of  interesting  new  stochastic  and  deterministic  mathe¬ 
matical  and  computational  methodologies. 
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We  also  note  that  the  approach  and  methodology  presented  here  are  useful  for  investiga¬ 
tion  of  a  wide  range  of  perturbations  other  than  disease  (e.g.,  loss  of  capacity  at  a  given  node 
such  as  a  factory  being  shut  down  for  some  reason)  in  supply  networks.  In  particular  they 
would  be  useful  in  the  assessment  of  risk  of  a  food-borne  pathogen  (e.g.,  salmonella,  listeria, 
etc.)  entering  the  food  chain  [27]  due  to  contamination  (either  accidental  or  deliberate)  at 
some  stage  of  the  supply  chain. 
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